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ABSTRACT 



^ ■ We investigate the motion of a particle around a low mass planet embedded 

in a non-turbulent gaseous disk. We take into account the effect of the gas 
structure that is modified by the gravitational interaction between the planet. 
We derive an analytic formula that describes the change of the semi-major axis 
of the particle due to the encounter with the planet using local approximation 
in distant encounter regime. Our final formula includes the effects of steady, 
axisymmetric radial gas flow, the global gas pressure gradient in the disk, planet 
gravity, and the structure of the gas flow modified by the planet's gravity. We 
compare the analytic results with numerical calculations, and indicate that our 
formula well describes the secular evolution of the dust particles' semi-major axes 
well, especially for small particles with large drag coefficient. We discuss the 
conditions for dust gap opening around a low mass planet and radial distribution 
of dust particles. Our formula may provide a useful tool for calculating radial 
distribution of particles in a disk around the planet. 

Subject headings: planet and satellites: formation — solar system: formation — 
celestial mechanics — methods: analytical 
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1. Introduction 

The motion of small particles in a gas disk has been investigated in detail since this 
is essential for planet formation theory. In a number of work, an axisymmetric gas disk is 
assumed and the effects of the structure of the gas disk on particle orbit are investigated, 
neglecting the effects of planet's gravitational force (Adachi et al. 1976, Weidenschilling 
1977). Their results show that particles with stopping timescale of the order of Kepler 
timescale migrate towards the central star rapidly, thereby imposing a serious barrier on the 
formation of planetesimals. 

Although this poses a serious question on the formation of planetesimals, there is another 
aspect of this fact in the formation of gas giant core. Gas giant core should reach the critical 
mass to accrete nearby gas rapidly (Mizuno et al. 1978), and its formation timescale seems 
to be very long (Pollack et al. 1993). If the global pressure gradient exerted on the gas 
disk will force the particles to migrate towards the central star, the particles at outer disk 
may reach the orbit of an already formed core and may accrete onto it. Kary et al. (1993) 
have shown that this process seems to be possible but most of the particles may miss the 
planet since they migrate very rapidly. Their discussion is based on the three-body problem 
including the effect of axisymmetric gas disk. 

The assumption of axisymmetric gas disk is a simplification in order to make the cal- 
culation tractable. The planet gravitationally interacts with the gas disk and changes the 
structure of the disk around it, producing the spiral density wave (see e.g., Goldreich and 
Tremaine 1979). Only recently does appear the investigation of particle motion around the 
planet embedded in a gas disk, which fully takes into account the structure around the planet 
(Paardekooper and Mellema 2004, 2006, Paardekooper 2007, Fouchet et al. 2007, Lyra et 
al. 2008). For instance, Paardekooper and Mellema (2004, 2006) and Paardekooper (2007) 
numerically calculated the particle motion around a high mass planet in the presence of a 
gas disk and discussed the accretion of the particles onto the planet and the formation of gap 
of particles of different size. They concluded that minimum mass of the planet that opens 
up a dust gap was 0.05M/, thereby indicating the possibility of detecting a low mass planet 
by future observations by observing a dust gap. 

Many of such works involve numerical studies on particle motion. In this paper, we 
present an analytic investigation of the motion of particles around a low mass planet em- 
bedded in a gas disk, taking into account of the structure around the planet. This is the 
full analytic solution of a particle motion around the planet in which the structure of the 
gas disk is fully taken into account. We consider a vertically averaged disk and perform 
two-dimensional analysis for simplicity. We use the results of (quasi-) linear calculation of 
the gas structure around the planet and calculate the distant encounter between a particle 
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and a planet embedded in a gaseous disk. We derive a formula which describes the secular 
change of the semi-major axis of the particle caused by the encounter with the planet using 
Hill's equations. It is not possible to describe the resonant capture discussed by Weiden- 
schilling and Davis (1985) or Kary et al. (1993); however, we take into account the structure 
of the gaseous disk that is modified by the planet, which is not considered in most of the 
previous analytic works. Our paper is complementary with previous numerical studies in two 
respects. Firstly, our analysis is focused on low mass planet, and secondly, our treatment 
is analytic. It is challenging to follow numerically the motion of dust particles embedded 
in a gas disk for long timescale (e.g., 10 5 Kepler time), solving gas and the motion of dust 
simultaneously. Therefore, analytic treatment may be necessary in order to understand the 
secular motion of dust particles. Moreover, basic physical processes in the problem become 
clear by analytic treatment. 

The plan of this paper is as follows. In Section [2j we show analytic treatment of the 
problem and derive the formula for the evolution of semi-major axis, equation (|68|) . This 
formula includes the effects of global pressure gradient of gas disk, steady, axisymmetric mass 
accretion flow of the gas, gravitational force from the planet, and the gas velocity structure 
modified by the planet. We then compare our results with numerical calculations in Sections 
|3] and HI and the limitations of our analytic approach are discussed. In Section we discuss 
the conditions for dust gap formation and long-term behavior of dust particles whose orbits 
are close to an embedded planet. Section [6] is for summary and future prospects. 

We expect that our results, equation (|68|) . can be used as a tool for estimate radial 
motion of the particle in a disk around the planet, since this is written in an analytic form. 
It may be possible to calculate the radial distribution of particles of various sizes using this 
formula and derive, for example, the opacity of the disk. We recommend readers who need 
only the final results of the particle motion just to refer equation (I68p directly and go to 
Section where we demonstrate the calculation of radial dust distribution around a planet. 

Before proceeding to the details of calculations, we note that there are two relevant 
length scales in the problem at hand: pressure scale height H of the gas disk and Hill's 
radius of the planet, th- Pressure scale height is an important quantity for the structure of 
gas modified by the gravitational force of the planet, while Hill's radius is important for the 
orbit of particles. Pressure scale height that is given by sound speed of the gas c divided by 
Keplerian angular velocity and Hill's radius of the planet are related by 



where r p is the semi- major axis of the planet, M p is the mass of the planet, and M c is the 
mass of the central star. For low mass planets considered in this paper, r-n/H is smaller than 
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unity. 

2. Analytic Consideration 

In this section, we analytically investigate the particle motion under the influence of 
gravitational force by the protoplanet and gas drag using local approximation. We show 
that particles with friction time of the order Kepler time will migrate towards the planet 
because of the gas drag, although gravitational scattering tends to repel the particle from the 
planet. This mechanism is different from the inward migration of dust particles caused by the 
velocity difference between the gas and dust [see e.g., Adachi et al. (1976) or Weidenschilling 
(1977)]. We also derive the effects of the modification of gas velocity due to the gravitational 
force of the planet. 

2.1. Basic Equations 

We use Hill's approximation to investigate the encounter between the particle and the 
planet. We consider the motion only in (x,y) plane where x is in radial direction and y is 
in the azimuthal direction. We assume the planet is fixed at the origin (x, y) = (0, 0) in this 
coordinate and neglect the mass of the particle. We write the velocity of the gas v g by 

3 

v g = — -Q p xe y + 5v g , (2) 

where Q p is the Keplerian angular velocity of the planet and e y is the unit vector in the 
y-direction. The first term represents the Keplerian motion around the central star and 
the second term Sv g represents the deviation from this motion, which, for instance, may be 
caused by the global pressure gradient of the disk. 

The equations of motion of the particle are 

x — 2Q p y = 3QpX — vx + F x (3) 

y + 2Q p x = -v(y+ hl p x^j + F y (4) 

where v stands for the drag coefficient. The quantity v is a reciprocal of the stopping time 
of a particle. Note that there are various notations for this parameter. In Weidenschilling 
(1977), stopping time is denoted by t e so v in this paper is equal to l/t e . In Adachi et al. 
(1976), stopping time corresponds to their A times gas density times relative velocity. We 
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write the sum of the force exerted by the planet and the part of the friction force caused by 
the non-Keplerian motion of the gas by F, 

F(x,y) = V-^£= + v6v g {x,y), (5) 

v x ~t~ y 

where G is the gravitational constant and M p denotes the mass of the planet. 

We give a brief note on the relationship between the dust size and drag coefficient. In 
this paper, we use the constant drag coefficient v for mathematical convenience. The dust 
size d and the drag coefficient v are related by v oc d~ x for Epstein law and v oc dr 2 for 
Stokes law. Simple relationship between the dust size and drag coefficient may be given by 



v 



1Q4 Po \ MV 1 J_ K 10-Vcm 3 



10" 9 g/cm / V lcm / lcm Po 

-2 _o 3 ^ 

in 1 1 d \ 10 V cm < d 



lcm J po 1cm 



where po is the background density and d is the size of the particle. Figure [T] shows the 
relationship between dust size and drag coefficient given by this equation for p = 10 _8 gcm~ 3 , 
10~ 10 gcm~ 3 , and 10~ 12 gcm -3 . 

The assumption of constant drag coefficient is not true for large bodies whose drag force 
is proportional to the square of the velocity difference. However, we expect that the results 
are qualitatively similar even for large bodies. 

We neglect the effects of gravitational potential caused by gas disk. It is a good approx- 
imation for non-self-gravitating disk as shown in Section 12.61 

In the absence of the planet and the friction force, the particle motion is the Keplerian 
motion around the central star. In our setup, this is represented by four orbital elements 
(b, h,k,y ), 

x(t) = b + r p h cos(f2 p t) + r p k sm(Q p t) (7) 
3 

y{t) — Vo — -^bVL p t — 2r p h sin(fi p t) + 2r p k cos(fi p t). (8) 

Physically, b is the difference of orbital semi-major axis between the planet and the particle, 
h and k are quantities related to the orbital eccentricity, and yo is related to the initial 
azimuthal position of the particle. 

We calculate the evolution of the orbital elements caused by the planet's gravitational 
force and the friction due to the gas. The time evolution of the osculating elements is given 
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by 

2 

b = ur v [h cos(fl p t) + k sin(ftpt)] + — F y (9) 

s Lp 

h = -vh l — [F x sin(fipt) + 2F y cos(f2 p i)] (10) 

k = -vk + — — [F x cos(fipt) - 2F y sin(jl p i)] . (11) 

Tpl Lp 



2.2. Change of Semi-Major Axis by Distant Encounter: General Treatment 

We now solve equations (}9~|)- (JTT]) . We assume that the particle is in circular orbit initially: 
h(t = — oo) = k(t = — oo) = 0. We derive the change of orbital semi-major axis Ab = b(t = 
oo) — b(t = — oo) by the encounter. 

The formal solutions of equations (fTUj) and ffTTl) are given by 

1 f°° 

h(t) = — / due~ vu {F x (t - u) sin[fi p (t - u)} + 2F y (t - u) cos[Sl p (t - it)]} (12) 

r p"p Jo 

i r°° 

k{t) = — — / due~ uu {F x {t - u) cos[fi p (t - u)\ - 2F y [t - u) sin[O p (t - u)}} , (13) 

r p"p Jo 

where F x (t) denotes F x (x(t),y(t)) where (x(t),y(t)) is the location of the particle at time 
t given by equations (j7]) and (jHJ) and so for F y . Substituting equations (fl2"j) and (ITBl into 
equation (jUJ), we find 

/oo 
6(t)dt 
-oo 

/-OO /-OO 

= — / dt / rfMe _I/M {F x {t - u) sin(fi p M) - 2F y {t - u) cos(j1 p w)} 

"p J-oo JO 

2 f 00 

+ — / dtF y (t) (14) 

^p oo 

Changing the integration variable t —>■ t — u in the first integral of the second line and using 
the formula 

due~ uu cos{Qpu) = / ra (15) 

V + sip 



and 



o 



d M e-^sin(fi pM ) = — f^— , (16) 
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we find 

This equation formally describes the amount of change of semi-major axis. We note 
that the integration with respect to t must be performed along the particle's trajectory. The 
force F consists of two parts [see equation (151) ]. One is the gravitational force of the planet 
and the other is the non-Keplerian motion of the gas velocity. We now see these effects 
separately. 



2.3. Planet Encounter 



First, we consider the change of the semi-major axis caused by the planet. We assume 
that the trajectory of the particle is close to circular orbit. We approximate F x by 

GM n 1 



F x sgn(6) 

With this approximation, we obtain 



b 2 (1 + (9/4)(y 2 ) 3 / 2 ' 
4GM n 



F x (t)dt = -sgn(b)- 



3b 2 



18) 



(19) 



For the integration of F y , we can use the well-known result of restricted three-body problem 
(see e.g., Goldreich & Tremaine 1980, Henon & Petit 1986, or Hasegawa & Nakazawa 1990). 
The result is 



F y (t)dt 



64 G 2 M 2 
243 b 5 tt 3 p 



1 2 



(20) 



where Kq and K\ are modified Bessel function of zeroth and first order. The derivation of 
this equation is outlined in Appendix [A] 



Substituting equations (1T91 and (1201 into equation (1T7I) . we obtain 



Ab = -sgn(6)4^ 



+ a- 



n 2 



b 2 V 2 + tt 2 b 5 V 2 + tt 2 



where rn is the Hill's radius of the planet defined by 



3M* 



1/3 



and a is a numerical factor 



a = 



128 



2K 



30.094 



(21) 



(22) 



(23) 
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The second term of equation ff2Tl) represents the effect of gravitational scattering. The 
particle and the planet tends to repel each other by mutual gravitational interaction, which 
is a well-known result of restricted three-body problem. This effect remains when drag force 
vanishes {y — ► 0). When there is no drag force, the gravitational interaction between the 
planet and the particle results in the excitation of eccentricity and the difference in semi- 
major axes of the two increases since Jacobi energy must be conserved in the absence of any 
dissipative force. 

When friction force is large, gravitational scattering is ineffective since the gas drag 
enforces the particle to move with the fluid element. The first term of the equation (1211) 
shows that the orbits of the planet and the particle tend to attract each other and this is 
most efficient when Q p ~ v. The intuitive explanation of this effect is as follows. When 
the particle feels the gravitational force of the planet, it is attracted towards the planet's 
position at first. Then, when the gas drag is effective, the gas enforces the particle to stay at 
the Kepler orbit at the location where the particle is attracted, resulting in the attraction of 
the semi-major axes of the particle and the planet. If the drag force is taken into account, 
the semi-major axis difference between the planet and the particle can shrink since Jacobi 
energy is not necessarily conserved. 

Since the attraction of the semi-major axes of the particle and the planet in the presence 
of drag force represented by the first term of equation (1211) is proportional to b~ 2 at large b, 
this overwhelms the effect of scattering represented by the second term, which is proportional 
to b~ 5 . Therefore, we conclude that particles far away from the planet are attracted towards 
the planet when the gas velocity equals Keplerian rotation velocity. 

In the vicinity of the planet, it is expected that gravitational scattering is more effective 
since this effect increases as b~ 5 . The value of b where these two effects are comparable is 
given by Ab = 0, that is 



This indicates that particles with v > f2 p may move towards the planet even in the absence 
of the global pressure gradient. 

2.4. Non-Keplerian Rotation of Gas Disk due to Pressure Gradient 

In this subsection, we consider the change of the semi-major axis caused by the effect of 
non- Keplerian rotation of the gas, which is due to the presence of global pressure gradient. 
The result of this section is already derived by previous studies such as Adachi et al. (1976) or 




(24) 
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Weidenschilling (1977), but we briefly show the results in this formulation for completeness. 

We parameterize the degree of non-Keplerian rotation by rj and write the velocity dif- 
ference 5v g by 

5v g = r]v p e y = const, (25) 

where v p is the rotation velocity of the planet around the central star, v p = r p Q p . The 
parameter rj is non-dimensional and its magnitude is the order of the square of the disk 
aspect ratio. This is negative for the disk with negative pressure gradient. 

Since 5v g is constant, it is easy to integrate (TTTj) to obtain 

Ab = 2 ^ T ^rk • (26) 

where T is the time taken for the particle to cross the box. As is well-known, this effect is 
most efficient for the particle with v ~ fl p . 



2.5. Steady Accretion Flow 

In this subsection, we consider the change of the semi-major axis caused by the steady 
accretion (or deccretion) flow. This may be modeled by the constant axisymmetric radial 
velocity 5v x . We parameterize the radial velocity by non-dimensional factor ( as follows, 

6v g = (v p e x = const. (27) 

The parameter ( denotes the ratio between the radial flow velocity and the Kepler velocity. 
If the dissipation of the gaseous nebula is due to the gas accretion, the sign of ( is negative 
and its magnitude is the order of 10 -6 , which is indicated from observation (see, e.g., Haisch 
et al. 2001). Drag force due to the gas is given by u5v g and therefore, using equation (TT7]) . 
we find 

A6 = ^ T ^- (28) 

We note that the steady accretion flow effect is important for small particles, which have 
large drag coefficient u, since such particles move in the same way as the gas flow. 



2.6. Non-Keplerian Rotation of Gas due to the Presence of the Planet 



In this section, we consider the change of the orbital semi-major axis of the particle 
caused by the deviation of the gas velocity from Keplerian rotation velocity due to the 
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presence of the planet. In this section, 5v g is caused by the gravitational interaction between 
gas and the embedded planet, in contrast to the previous sections. We show that only 
axisymmetric structure of the disk causes the change of the semi-major axis of the particle. 

In order to calculate the velocity perturbation caused by the planet, we solve vertically 
averaged Euler equations using local shearing-sheet approximation: 

— + V-(S«) = 0, (29) 
dv c? 

— + (v • V)« = 3Q 2 p xe x - -VE - V^ p - 2Q p (e z x v), (30) 

where £ denotes the surface density, v is the gas velocity, c is the sound speed, and f2 p is the 
Keplerian angular velocity of the planet. We assume an isothermal disk, where c is constant, 
for simplicity. The gravitational potential of the planet, ip p is given by 

^ = / 2 ^ 2 ^ 2 ' ( 31 ) 

a/ x l + y z + e z 

where e is the softening parameter. Taking the rotation of the equation of motion, we obtain 
the equation for vorticity, 



9 



1 (dVy 8V X | ^ 

S \ dx dy r 



0. (32) 



For small mass planets we are interested in this paper, it is possible to calculate the 
velocity perturbation by linear analysis. However, as shown later, it is necessary to calculate 
the flow up to the second order in order to obtain the correct results for the particle motion. 
We assume that background surface density S is constant and background gas is rotating 
with Kepler velocity v = —(3/2)Q p xe y . We neglect the effect of the global pressure gradient 
or steady mass accretion since it is calculated separately in the previous subsections. 

The first order perturbation of the flow is caused by the planet gravity. We denote 
the first order perturbation of physical variables by superscript (1). The linearization of 
equations (1251) and (13U1) are 
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Assuming that vorticity is not created by the formation of the planet, we have 

d „ m d „ m 1 5E« 



0. 



(36) 



see, e.g., Narayan et al. (1987). Note that E appears only in the form of SYr- 1 ' /E and 
therefore, this value itself does not depend on the background surface density. The source 
term of these equations is given by planet potential, ipp so the perturbation is proportional 
to the planet mass M p . For homogeneous (ip p = 0) equations, only dimensional parameters 
are c and Q p . Taking fi" 1 as a unit of time and H = c/Q p as a unit of length, homogeneous 
equations become scale free. Since the amplitude of the perturbation is proportional to the 
source term, which is ip p in this set of equations, the perturbation amplitude scales with 
ipp/c 2 , in other words, 

<5E« GM n 



oc 



E -■ He* (37 > 
as shown by Tanaka et al. (2002). Once we know the perturbation amplitude of a specific 
value of GM P / He 2 , we can easily obtain the actual amplitude of density perturbation 5E for 
different background density, protoplanet mass, sound speed and so on. Linear calculation 
shows that the coefficient of proportionality of equation ( |371) is the order of unity. 

The first order perturbation becomes a source of second order perturbations through 
non-linear terms of basic equations. We denote the second order perturbation by superscript 
(2). The set of second order perturbation equations is 



d_ 

Iff 



3 A 

2 pX dy 
d 



ox V 



+ 7T' Sl 
ox 



d 



(2) 



+ 



d_ 

dy 



Sv 



(2) 



+ 



dy 



(3* 



dt 2 pX dy 



SEW d SEW 

En OX En 



2%8v { y ] + c 2 



ox 







L-lo JL 
dt 2 pX dy 



5v ( y 2) + ^-Q p Sv^ + c 2 



d SEW 
dx E 

dy 
d (5EW 



(39) 



9 5E« d 8YF> 

= c 

E dy E 

The second-order vorticity equation is 

dy Vx dx° y + 2 iZp E S 



dy S 



dy 



JEW 
S 



(40) 



(41) 
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Since the first order perturbation scales with GM P /Hc 2 (see equation fl371) ). it is clear that 
the second order perturbation scales with (GM p /Hc 2 ) 2 . We assume that the perturbation is 
stationary with respect to the planet so in the corotating frame with the planet considered 
in this subsection, d/dt = 0. 

We now see equations (!33|) - (l35|) in more detail to find an analytic expression of Ab. We 
first give an order-of-magnitude discussion to show that the gravitational force produced by 
the density fluctuation is small compared to the gravitational force by the planet. We expect 
that the order of magnitude of the gravitational force by the spiral density fluctuation is 

G5ZH 2 

Spiral ~ — jp— , (42) 

since the characteristic length scale of the spiral is of the order of the scale height. On the 
other hand, the gravitational force by the planet exerted on the particle at the distance of 
the order of the scale height is 

GM p 

-fplanet ~ ~Jp~' ^3) 

Using equation fl37j) . we expect that these two forces are related by 

#^ ~ ^ - Q'\ (44) 

^planet Ciip 

where Q is Toomre's Q parameter of the disk. Since protoplanetary disks in planet forming 
phase are expected to be gravitationally stable in general, Q ^> 1. Therefore, we conclude 
that the gravitational interaction between the particle and the disk is negligible. We later 
confirm this numerically. In this section, we consider only the effect of the velocity fluctuation 
of the gas disk in the presence of the planet. 



2.6.1. First- Order Axisymmetric Mode 

We first consider axisymmetric modes, or azimuthally averaged quantity. For any per- 
turbation quantities 5f, we denote azimuthal average by bars, 



1 



Ly/2 



Sf(x) = — Sf(x,y)dy. (45) 

L y J —Ly/2 

Using Euler equations ([55i) - ([551) and vorticity equation flSSD, we can derive 



SvP = 0, (46) 



<5E« 2 d 



dx 2 6V V C 2 6V V ~ 2c 2 dx & 
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For x <C L y , azimuthal average of dip/dx is given by 



dip 
dx 



1 



GM p x 



y J-oo [X 2 + y 



2^3/2 



dy 



2GM P 

LyX 



(49) 



where the integration range is extended to infinity. The appropriate boudary condition is 
such that the perturbation vanishes for \x\ — > oo. Since homogeneous solution of equation 
HJ) is exp[±x/H], Green's function of the differential operater of equation (1481) is 



G{x]x') = - — 



-(x-x>)/H d ( x _ ^ + e (x- x >)/H d ( x , _ x j 



(50) 



where 9(x) is step function. Therefore, the solution of equation (148]) is given by 



8v 



(i) 



H 2 tt p GM P 
2L y He 2 



e -(*/H) E i 



- e x/H Ei - 



x 
H 



(51) 



Using equation (T47j), the surface density perturbation is given by 



5EW 



H GM P 

Ty HC 2 



e-(^)Eif-) +e x ' H Ei(-- 



H 



x 



H 



(52) 



where Ei(x) denotes exponential integral. Equations (!5"T!) and (152|1 are valid for locations 

x <ti L y . Figure [2] shows the profile of L y 8vy\x)/c normalized by (GM P /Hc 2 ) and compares 
this with analytic expression (j5ip . Difference between numerical and analytic results is about 
25% at x ~ 5if. We have checked that this error significantly decreases when box size in 
x-direction is smaller than that of y-axis. 



2.6.2. Second- Order Axisymmetric Mode and Path Line of the Flow 

The second-order axisymmetric modes are obtained by equations (139]) . fj40|) . and (|4ip . 
We note that equation of continuity fl38l) and the ^-component of equation of motion (j40l is 

not independent for axisymmetric modes. We use only dVx later. This is given by 

W=~Sv^Svi 1) . (53) 

We then consider the path line of the flow. This gives the motion of the fluid element 
and therefore gives the motion of the particle with v — > oo. Since we are interested in the 



-14- 



change of the orbital semi-major axis of the particle, we consider the motion of the fluid 
element in the x-direction. Path line is given by 



dx 

— = 5v x (x(t),y(t)) (54) 



dy 
dt 

where v c is the unperturbed velocity 



v c + Sv y (x(t),y(t)) (55) 



3 

v c = —topX (56) 

and 8v x and 5v y include both first and second order perturbations; for instance, 5v x = 
Svi 1] + 8vi 2) . 

The path line of the unperturbed flow is given by 

x c {t) = b = const (57) 
Vcit) = h- 3 - b n p t (58) 

where we assume that the fluid element is at (b, L y /2) at t — 0. This is simply a Keplerian 
circular motion. We denote the time when the particle reaches the other end of the box in 
the unperturbed flow by t = T, where 

T = (3/2)MV (59) 

We solve equations (154")) and (|55|) using perturbation methods and obtain the change of 
x-coordinate, Ax, of the fluid element after it crosses the box 

Ax = x(T) - x(0) = f dt6v x (x(t),y(t)). (60) 

Jo 

We divide the motion of the fluid element into the unperturbed motion and perturbation 

x{t) =x c + 5x{t) (61) 
y(t)=y c (t) + 5y(t). (62) 

Expanding right hand side of equation (1601) upto the first order of 5x and 5y and using 
equations (I54T) and (|55l) . we obtain 

Ax= / 5vf\x c ,y c {t))dt+ dh-^-ix^y^)) dt 2 5v£\x c ,y c {t)) 

r T 8Sv (1) r h 

+ dti^-(x e ,y c (ti)) dt 2 6vV>(x e ,y e (t)) ( 63 ) 

Jo oy Jo 
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where we have used 8v x ^ = and take the terms up to the second order. Using the first order 
vorticity equation (|36|) and equation of continuity, it is possible to show that the right hand 
side of equation (I63I) vanishes, see Appendix [B] Therefore, we conclude that fluid element 
returns to the original radial position after crossing the spiral density wave, 

Ax = 0. (64) 

We note that second-order perturbation is essential to derive this conclusion and therefore, it 
is necessary to consider the fluid motion upto the second order in finding the correct motion 
of the particles in a disk. 

The fact that the fluid element does not move in the radial direction indicates that there 
is no radial mass flux up to second order perturbation. It is possible to show this directly 
by calculating Hv x , 

J dyEv x = J dyE 8v^ + S J dy^—Sv® 

= ~^ So / w^m 6 ^ + So / dy5 ^ 6v ^ 
= ~W< W1))2 

= 0, (65) 

where we have used equation (|53|) in the second line, vorticity equation (|36|) in the third 
line, and periodic boundary condition in the last equality Q. This is in contrast to the case 
of the sound wave propagating in a static, homogeneous medium, which carries the linear 
momentum and therefore mass flux is present (Landau and Lifshitz 1959). We demonstrate 
this in Appendix ICl 



2.6.3. Ab Caused by Spiral Density Wave 

We are now in the position to calculate Ab due to the structure of the gas disk modified 
by the planet's gravity. Using equation (1171) and drag force is given by u8v g , we have 

v 2 r°° IvVL f 00 

Ab = 7^Tm Sv A<t),y(t))dt+-—f- Sv y (x{t),y(t))dt. (66) 



1 It is actually possible to show that, in the shearing-sheet approximation, the radial mass flux vanishes 
for all orders of perturbation in the same manner as presented here. 
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The integration must be performed along the path of the particle. The first term of equation 
f l66|) dominates when v — > oo. In the perfect coupling limit, since the particle traces the flow 
of the fluid element, this integral should vanish. We expect that this integral always remains 
small and set the first term of equation f[6T)j) to be zero. We expect that this approximation 
is valid since the second term that is proportional to uQ p / [y 2 + Q 2 ) dominates Ab when v 
is small. We later confirm this numerically. 

We then consider the second term of equation (1661) . The leading order of this integration 

is obtained by approximating that the particle orbit is circular. Since Sv y 1} ^ 0, it is enough 
to consider this contribution. Therefore, we arrive at the following expression of the change 
of orbital semi-major axis caused by the spiral structure of the gas 

4 1 vVL 



where Sv y is given by equation (1511) . 

We note that in deriving equation (l6Tj) . we have repeatedly used the assumption of no 
vorticity source, equation (13"6"j) . In the presence of vorticity source, 8v y may be different 
and the mass flux may be present. 



2.7. Analytic Expression of the Change of the Orbital Semi-Major Axis 



Adding equations fl2T|) . (126]) . fl28l) . and fIBTl) . we obtain the expression for the secular 
evolution of the orbital semi-major axis of the particle. We write in the form of the rate of 
the change by dividing Ab by T = L y /(3/2)Q p \b\, where L y is the box size of the y-direction 
of the coordinate system considered. Time T may be interpreted as the time interval between 
successive conjunctions between the particle and the planet. The result is 
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(68) 



where rn is Hill's radius of the planet, the definitions of r\ and ( are given by equations ([25 
and (12T|) respectively, and the value of a is 30.094. 



If the planet mass is sufficiently small, the leading contribution comes from the first 
term, which is due to the global pressure gradient exerted on the gas disk. Comparing the 
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order of magnitude of the first and the third term of equation fl68|) . the third term becomes 
dominant for particles with 



where we have assumed L y ~ r p , rj ~ (H/r p ) 2 and all the numerical coefficients are neglected. 
For low mass planets which does not form a gap in the gas disk, Hill's radius is generally 
smaller than the scale height. Therefore, the motion of the particle is predominated by the 
effect of global pressure gradient. However, for large planets with mass of order Jupiter, 
the effect of the global pressure gradient is small compared to the gravitational interaction 
between the planet, and equation (I68I) indicates that if the profile of the gas disk is neglected, 
the particles may be attracted towards the planet. This seems to be consistent with what 
is found by Paardekooper (2007). He found that when axisymmetric gas disk profile is 
considered, plants with the mass of the order of Jupiter would accrete the particles with 
stopping time comparable to Kepler time, while accretion of such particles around planets 
with 0.1 times Jupiter mass would be inefficient (see figures 7 and 11 of Paardekooper 2007). 
When global pressure gradient of gas disk drives the particles to migrate towards the central 
star, a part of particles with v ~ Q p may pass the orbit of the planet without accumulating 
onto it, since the particle migration rate is very fast and the perturbation of the orbital 
elements of the particles is not related to the planet location (Kary et al. 1993). If the third 
term of equation (|68|) dominates over the global pressure effect, on the other hand, particles 
will accumulate onto the planet since their orbital elements are strongly perturbed at the 
location of the planet. 

If the disk model with flat pressure profile is considered, the behavior of the particle 
around the planet is complicated. Figure [3] shows the contour plot of equation (168]) for the 
planet mass GM p /Hc 2 = 10~ 2 , flat pressure profile (rj = 0), and no mass accretion onto the 
central star (£ = 0). For small particles [y — > oo), there is no secular motion of particles since 
all the terms of equation (168|) vanishes if ( = 0. The motion of small particles is determined 
by the strength of accretion of the background flow. For large particles (v — > 0), on the other 
hand, the fourth term, which is the gravitational scattering by the planet, is dominant and 
the particles are always repeled away from the planet. For particles with intermediate size 
(v ~ Jlp), all the terms must be considered. Particles located sufficiently far away from the 
planet are scattered away because of the contribution from the last term of equation (168]) . 
which is due to the axisymmetric mode of the perturbed flow structure caused by the planet 
gravity. Particles at the intermediate distance (b/H ~ 1), on the other hand, are attracted 
towards the planet because of the third term. If the particles are close to the planet, they 
are again scattered away because of the contribution from gravitational scattering. 




(69) 



The equilibrium distance between the third and fourth term is given by equation (124|) . 



-18- 



and this becomes close to the planet with increasing v, and it eventually becomes smaller 
than the Hill's radius, where we expect that the particle motion is largely dominated by 
the planet's gravity and our treatment of distant encounter breaks down. Therefore, we 
conclude that the particles with intermediate but still smaller drag coefficient, on one hand, 
will accumulate around the planet, but not fully accrete onto the planet. On the other hand, 
particles with larger drag coefficient will accrete onto the planet since gravitational scattering 
is ineffective in this case. Equating the equilibrium distance given by equation (1211) to Hill's 
radius, the critical drag coefficient separating these two regimes is of order 10O p . 

We note that since equation fl68|) is obtained under the assumption of distant encounter, 
it does not accurately describe the dynamics of particles located within the Hill's sphere of the 
planet. In Section [31 we perform a numerical calculation of Hill's equations and investigate 
the validity of assumptions in the derivation of equation (|68j) . It needs to be checked whether 
particles will really accrete onto the planet when they approach towards the planet. For close 
encounter, the detailed structure of the gas around the planet must be considered. Inaba & 
Ikoma (2003) have shown that the capture cross section of particles by the planet may be 
enhanced when the atmospheric structure is considered. 



3. Numerical Calculation in Local Coordinate 

3.1. Numerical Methods 

To investigate the accuracy of the description by equation ( 1681) for the evolution of 
orbital element of a particle whose semi-major axis is close to the planet embedded in a 
gaseous disk, we perform a series of numerical calculations and compare the results with 
equation (1681) . 

In our calculation, we solve Hill's equations ([3]) and (0J using fifth-order Runge-Kutta 
method with variable time step. We are primarily interested in the effects of spiral density 
wave and the gravitational force by the planet centered at the coordinate system: we neglect 
the effect of the global pressure gradient and mass accretion, r\ = ( = 0. The velocity field 
of the gas is therefore 

5v g = 5v {l) + 5v i2 \ (70) 

where 5v^ and Sv^ are first and second order perturbations, respectively, described in 
section [231 We also note that for gravitational force exerted by the planet, we use 



GM P 

V x2 + y 2 



(71) 



19 



We do not use the time-averaged force given by equations (IA8I) and (IA9I) . 

The modification of gas velocity due to the gravitational force by the planet is obtained 
by second-order perturbation analysis. We numerically integrate equations (153|) - fl5o'j) for the 
first order solution and (138]) - (140]) for the second order solutions. We make use of Fourier trans- 
form methods presented by Goodman and Rafikov (2001). This method is particularly useful 
in calculating the gravitational potential produced by the density fluctuation. We fix the 
spatial resolution (therefore, the maximum wave number in Fourier space) (Ax/H, Ay/H) 
to be (9.8 x 10 -3 , 7.8 x 10~ 2 ), which is enough to resolve the effective Lindblad resonance (see 
e.g., Artymowicz 1993). We fix the box size in the ^/-direction such that —80H < y < 80H 
and impose a periodic boundary condition. For the x-direction, the resolution in Fourier 
space in k x direction is varied from 10 4 to 10 6 according to k y in order to resolve very fast 
oscillation at k x ^> k y . We also use the linear window function given by Goodman and 
Rafikov (2001). We set softening parameter e to be 0.01H. 

The box size of our calculation in ^-direction is —80H < y < 80H, to be consistent 
with the calculation of the modification of the velocity field. In integrating Hill's equations, 
we impose a periodic boundary condition in the y-direction and integrate the equations of 
motion until t = SOOOfip 1 . By periodic boundary condition, we mean that when the particle 
has reached the boundary of the box in y-direction, y = ±80H, it is reintroduced into the 
other end of the box, y = ^80H, while the values of x, v x , and v y are kept unchanged. 
In principle, Hill's equations describe the motion of the particle only in the vicinity of the 
planet, and hence there is no reason that the phase of the epicyclic motion is conserved when 
it is reintroduced. However, we expect that if the box size is the same as 27rr p , which is the 
circumference of the disk at the orbit of the planet, this periodic boundary treatment may 
simulate qualitative properties of the orbit in the global calculations. The condition 



is realized, in the parameter of our calculations, if we consider a disk with the aspect ratio 
H/r p = 2vr/160 ~ 0.04. 

Once we know the data of (x, x, y, y), we can calculate the osculating elements (b, h, k) 
by, using equations (J7|) and (jSJ), 



L y = 27rr p 



(72) 



b = Ax + —y, 



(73) 
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(74) 




(75) 
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In the following sections showing results of calculations, we refer to (b, h, k) using these 
relations. These values are not averaged over the encounter, but are defined at any time of 
the orbit. 

The initial condition of the calculation is that the particles are located at the edge of 
the box and moving in a circular motion. When a particle approaches to the planet and 
the mutual distance becomes smaller than half of the Hill's radius of the planet, we stop 
the calculation since we are primarily interested in the distant encounter. Such particles 
should be trapped by the planet's gravity and eventually accrete onto the planet. In order 
to obtain a smooth results in the runs with large drag coefficient, we make use of smoothing 
methods described in Appendix [D] in calculating the perturbed gas velocity at the location 
of the particle from the data obtained by hydrodynamic equations. We fix the background 
density to be 10~ 10 g/cm 3 , which is an order of magnitude smaller than the Minimum Mass 
Solar Nebula (Hayashi et al. 1985). Toomre's Q parameter of the disk at 1AU with is 

Q- 1 ~ ^ ~ 10- 4 (76) 

We assume the mass of the planet as GM P /Hc 2 = 10~ 2 that corresponds to 0.2M e for 
c = 10 5 cm/s and H = 0.05AU. In the following sections, we show the results without the 
gravitational force by gas. We have checked that gas gravity does not affect the results in 
our parameter range. We also neglect the effect of global pressure gradient and steady mass 
accretion in our calculation, r\ = ( = 0, in order to see the effect of spiral wave and planet's 
gravity. 



3.2. Results 

3.2.1. Properties of Orbital Evolution 

We first review the properties of orbital evolution of a particle with zero friction force, 
which is obtained by setting v = 0. Figure H] shows the evolution of b for a particle initially 
located at b = H. The values of b are obtained according to equations (1731) when the particle 
is at the box boundary. Also indicated in Figure H] is the variation of b for < t < 250f2~ 1 . 
The values of b varies rapidly when the particle encounters with the planet, when y ~ 0. 
We note that since the box size is 160H, time taken for the particle to cross the box is 
~ 160/(3/2) ~ 100 Kepler times. 

It is indicated that the semi-major axis difference b first increases but then shows an 
oscillation with the period ~ 1200 Kepler times. This oscillation is caused by the excitation 
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of eccentricity by the perturbation due to the planet. If the particle comes into the box with 
finite eccentricity, the distance between the particle and the planet at the closest approach 
is different depending on the phase of the epicyclic motion, and the semi-major axis of the 
particle can increase or decrease. The period and amplitude of this oscillation depend on 
the location of the particle. If the phase of the epicyclic motion among successive encounters 
between the planet and the particle matches, the period of the oscillation becomes longer 
and the amplitude becomes larger. The condition for strong amplification is given by 



Figure El compares the evolution of semi-major axis of the particle starting at I = 5 position 
(b = 3AH) and at b = 3.8H. It is clear that the particle at I = 5 position experiences the 
long-term, strong oscillation of semi-major axis. 

Equation (I77j) indicates that there is a resonance between the box crossing time and 
the epicyclic motion of the particle, and this appears as a result of our treatment of periodic 
boundary in Hill's coordinate. Thus, in general, we should not consider this resonance to be 
physical if we arbitrary choose the box size in the y-direction. However, if we take the box 
size L y as the circumference of the disk, 27rr p , the crossing time T corresponds to the period 
of synodic encounter. Therefore, it may be possible to interpret this resonance as + 1) 
resonances in global problems. In section HJ we explore the correspondence between the 
Hill's approximation and the full global problem. 

We now consider the problem with gas drag. Figure El shows the evolution of b for 
particles with various drag coefficients initially located at b/H = 1 as a function of time. 
The values of b when the particles reach at the edge of the box are plotted. Particles with 
drag coefficient larger than u/Q p = 10~ 2 shows a systematic change of semi- major axis while 
that with v/Q p = 10~ 4 shows a systematic change plus oscillation. In a simple restricted 
three-body problem, where drag coefficient v equals zero, semi-major axis does not show a 
systematic change but only oscillation, as discussed in previous paragraphs (see also Figures 
IU and The gas drag causes a systematic change of Jacobi energy and thereby results 
in a systematic change of particle's angular momentum. The timescale during which the 
oscillation of semi-major axis damps is given by ~ v^ 1 . Since the oscillation of the semi- 
major axis comes from the excitation of eccentricity during the synodic encounter, this 
oscillation does not appear when the eccentricity is damped during successive encounters. If 
the eccentricity is damped during successive encounters, the assumption of initially circular 



fipT = 2-kI, 



(77) 



where I is an integer and T is the time taken to cross the box in y-direction 



T 



(3/2)6fip" 
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orbit in deriving equation (168]) is valid for every encounter. Therefore, although we have 
computed only one encounter, it is possible to use equation (1681) to predict the long-term 
evolution of the orbital semi-major axis especially for small bodies with large drag coefficient. 



3.2.2. Comparison between Analytic and Numerical Calculations 

We now show to what extent equation (1681) describes the evolution of the orbit of a 
particle encountering with a planet embedded in a gaseous nebula. Since the particles show 
a systematic change of semi-major axis, we take the average radial velocity of the particle 
by calculating 

&(*fin) ~ &(*ini) ( 7Q ) 
^calc 

where t- m \ = is the initial time, t^ n is the last time when the particle reaches the box 
boundary before the calculation is stopped at t — 5000fi _1 , and t ca i c = t^ n — t^- The 
precise value of t^ n is obtained by numerical calculation, and it varies with the value of b. 
Semi-major axis of the particle at each time is calculated by equation (1751) . 

Figures [7J19J show the results of long-term evolution of semi-major axis of particles ini- 
tially located at various position with u/Q p = 10~ 2 , 1, and 10 2 respectively. Horizontal axes 
of Figures [7J19J are the initial semi-major axis of the particle, b(ti n {), and the vertical axes 
show the average rate of change of semi-major axes obtained by equation (1791) . We also plot 
our analytic formula given by (|68|) . It is clear that the analytic formula (|68p describes the 
results of numerical calculations well, at least quantitatively, for all values of drag coefficient. 
Deviation from the analytic value for u/Q p = 1 at b/H ^ 2 comes from the deviation of 

the data of 5vy from the analytic value shown in Figure [2] while for small 6, we expect 
that the numerical results deviate from the analytic calculation since eccentricity effect is 
not negligible. The match between numerical calculation and analytic formula is better for 
particles with large drag coefficient, since their eccentricities remain smaller. We think that 
the deviation from the analytic value for v/Q p = 10 2 at b/H ^ 3 is due to smoothing. 

One significant difference between analytic results and numerical calculation is the spike- 
like structures that appear in the calculations with v/VL p = 1CT 2 . This structure is due to 
the resonance effect given by equation (1771) . We show in Figure [7J the location of resonance 
with I = 4, 5, and 6. It is clear that Ab/T shows a spike at this location. We have not 
found such resonance structure for calculations with v /Q p ^ 1. This is because eccentricity 
is sufficiently damped between successive encounters. 

Another interesting result is the motion of the particle with Q p ~ v. Figure [10] shows 
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the evolution of the semi-major axes of the particles with u/Q p — 1. It is clear that the 
particles migrate towards the planet and stagnate at b ~ 0.35H. Since we set rj = ( = in our 
calculation, the orbital change of the particle is caused by the third, fourth, and fifth terms of 
equation (|68p . The last term is negligible for small b compared to the third term. Equation 
(|68|) indicates that these particles accumulate near the planet's orbit where gravitational 
scattering represented by the fourth term of equation (1681) balances with the attraction 
towards the planet represented by the third term, which is the result of combination of the 
planet's gravity and the gas drag, as discussed in Section |2TB"1 The condition for the balance 
between these two effects is given by equation (liMj) . For chosen parameters in numerical 
calculations, rn = 0.15H so the equilibrium distance is given by ~ 0.3H, which is close to 
the value obtained in numerical calculation ~ 0.35H. We consider that the small difference 
comes from the effect of eccentricity and close encounter. This effect seems analogous to 
"the shepherding effect" . We note that this effect appears when the first term of equation 
(|68|) . Ab by global pressure gradient, is small compared to the third term. This condition is 
expressed in terms of 77: 

N<^4~^- (so) 

1 " b rl b M* v ; 

If rj is of the order of square of disk aspect ratio, as in standard parameters, this condition 
is rewritten as 

For Mp/M* = 1CT 6 and H/r p = 0.05, this gives b/r-^ <^ 0.2, which is too close to the planet 
and the assumption of distant encounter is violated. Kary et al. (1993) has shown that, in 
the presence of global pressure gradient, a particle with v ~ Q p falls onto the central star so 
rapidly that it bypasses the planet without being trapped. 



3.3. Limitation of Analytic Formula 

Equation (1681) is derived by assuming that the deflection angle of the particle by each 
encounter with the planet is small (distant encounter) and that the particle returns to a 
circular orbit before the next encounter. 

First, the assumption of initially circular orbit breaks down if we consider the particle 
with small friction. The oscillation of semi-major axis observed in figure [6] for particles 
with v = 10~ 4 f2 p indicates that the damping of eccentricity between successive encounters 
is insignificant for these particles and that the finite eccentricity effect should be taken into 
account. Quantitatively, the period of synodic encounter must be longer than the particle's 
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stopping time. Therefore, our treatment is limited in the cases with 



""V X -, (82) 



(3/2)6fi p ~ u 

where the left hand side is the synodic period and the right hand side is stopping time. 
Rearranging this condition, we have 

b Arc r n v , . 

£ ^~E7T- (83) 



H 3 H Q p 



Adopting standard parameters, we have 



A^of^iVV^, (84) 



H ^ \H/r } 

We now consider the assumption of distant encounter. If the particle's semi-major axis 
is too close to the planet, the assumption of small deflection angle, or distant encounter, 
breaks down. In Figure [TTJ we plot the amount of deflection of the particle's semi-major 
axis after the very first encounter with the planet. We quantify the amount of deflection by 

1*01 ' (85) 

where t\ = = is the initial time, and £2 ~ L y /(3/2)bQ p is the time when the particle 
reaches the box boundary after the first encounter. The precise value of t 2 is obtained by 
numerical calculations. Particles in black regions are trapped by the planet, that is, in our 
criterion, mutual distance between the planet and the particle becomes smaller than half of 
the Hill's radius, or in horseshoe orbit. Note the similarity of the shape of the contour of 1% 
deflection with the contour shown in Figure [3J 

In three-body problem with zero drag force, it is known that the assumption of small 
encounter is violated for particles with b <^ 2.5rn if we consider the very first encounter. 
Our calculation is consistent with this for particles with small drag, v <C Q p . However, for 
particles with large drag, v ^> Q p , it can be seen that this condition is relaxed and it is 
indicated that the assumption of small deflection angle can be used even down to b ~ th- 
Figure [T2l shows the orbits of particles initially located at b = 2r^. Orbits of particles with 
v/Qp = 1CT 4 and 10 3 are shown. Particles with very large drag coefficient can remain in 
the orbit close to the initial semi-major axis since gas drag prevents it from falling onto the 
planet. 

After multiple encounters, the condition of distant encounter is violated for particles 
with b <^ 2.5 — 2y / 3^H i n three-body problem without gas drag. The contour of 1% deflection 
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in Figure [TT1 passes close to 2a/3Vh i n the limit of small drag. Therefore, we infer that this line 
gives the limitation of our analytic treatment assuming small deflection. This assumption is 
better for particles in wider range of semi-major axes for larger drag. 

4. Comparison between Local and Global Calculation 

So far we have performed analytic and numerical calculations of the orbits of particles 
using local coordinate systems. It remains to be checked how our local treatment simulates 
the real orbit in global coordinate. 

We have performed a global calculation with the gas rotating at Keplerian rotation 
velocity where we have neglected the effect of the spiral pattern around the planet, global 
pressure gradient, and gas accretion towards the central star, rj = ( = 0. We solved the 
equations of motion in a frame rotating with the planet and integrate until the particle 
encounters with the planet for 100 times. We calculated the evolution of the semi-major axis 
of a particle at various initial semi-major axes and obtain Ab/T over the calculation time. 
In the following, we show calculations in which the particle is located at the opposition point 
initially, but we have checked that the initial locations do not affect the results as long as we 
consider particles with u/Q p > 10 -2 . Particles that reach the distance with the planet less 
than Hill's radius are removed. We set disk aspect ratio H/r p = 0.04, which gives 2nr p = L y , 
where L y = 160H is the box size we have used in calculations with Hill's equation. For planet 
mass, we assumed M p /M* = 6 x 10~ 7 , which corresponds to GM P /Hc 2 = 10~ 2 , the value we 
have adopted in calculations with Hill's equations. In order to compare the results of global 
calculations with Hill's equations, we also integrate Hill's equations ([2]) and fl3| by setting 
5v g = 0. 

Figures [TBlfTBI compare the results of Ab/T obtained by global and local calculations. 
Horizontal axes of Figures [1311151 show the initial semi-major axes of the particle and the 
vertical axes show the average rate of change calculated by equation (1751) . In global calcula- 
tions, 6(tfi n ) is the semi-major axis difference between the particle and the particle calculated 
when the particle has reached the opposition point after 100 encounters. 

For all parameters, global and local calculations show a reasonably good agreement. 
Relative difference between global and local runs is most significant for particles with u/Q p = 
1. We infer that this is because in the case of v/VL p = 1, the radial motion of the particle 
is the fastest, and the effect of curvature is most effective. In runs with u/Q p = 10 2 and 
u/Q p = 10~ 2 , radial motion of the particle is smaller by one or two orders of magnitude than 
calculations with u/Q p = 1. 
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In runs with v/Q p = 10~ 2 , we observe a disagreement in resonance locations, since 
curvature effects are not taken into account. However, since the disagreement is less than 
10%, it may be still possible to say, for example, that I = 4 resonance in Hill's coordinate 
mimics 4 : 5 resonance in global coordinate. Note that the relative strength of the resonance 
matches well with global runs. 

Global analysis of j/ (j + 1) resonances including drag force is done by Greenberg (1978). 
It is shown that eccentricity and the angle between the longitude of conjunction and the 
longitude of pericenter of the particle converges to a fixed value with the timescale of z/ -1 , 
but depends on the particle's semi-major axis. The value of eccentricity is maximum at the 
resonance. We have observed similar behavior in local calculations. It is actually possible 
to do analytic calculation similar to Greenberg (1978) in Hill's system and derive orbital 
evolution including the effect of resonance. We do not show it here since it is necessary to 
do all the analyses presented in Section [2] again including the effect of finite box size. We 
simply note here that this resonance effect does not appear in our formulation presented in 
Section [2] since we have neglected the effect of finite length in the ^-direction in deriving 
equation (1681 . 

In runs with u/Q p = 1, particles initially at b <^ 1AH shows the decrease of the average 
rate of semi-major axis decay. This is because the particle is trapped in the orbit in the 
vicinity of the planet as a result of shepherding effect discussed for calculations with Hill's 
equation in Section 1X2"! Figure [TBI shows the evolution of the semi- major axis of the particle. 
Osculating elements at the opposition points are plotted. Comparing with Figure [TD1 it is 
clear that the shepherding effect discussed in local calculations also exists in global runs. We 
note that the amplitude of oscillation after particles are trapped is larger in global runs than 
local calculations, and the location of the shepherding orbit is slightly closer to the planet 
than the local calculations. Note in passing that in equation (4.11) of Adachi et al. (1976), 
it is indicated that for r] = 0, damping of eccentricity causes a decrease in semi-major axis. 
However, Adachi et al. (1976) do not take into account the effects of the planet's gravity. 
The shepherding presented here is entirely caused by the presence of the planet. 

By comparing global calculations with local runs, it is indicated that as long as the 
orbits close enough to the planet are considered, local calculations reproduce the results of 
global calculations well even if curvature effect is neglected. However, the effect of constant 
box size in the ^-direction causes the small difference between the global and local runs. This 
is indicated by the mismatch between the locations of the resonances. It is also indicated 
that the local calculations reproduce the results of global calculations most efficiently for 
small particles, such as v/Q p ^ 10 _1 — 1CT 2 . Local approach that is used in Section[2]is valid 
when semi-major axis difference b is much smaller than the scale of the semi- major axis of 
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the planet: 

b , , 

— < 1, (86) 

which is equivalent to 

7? « T f m 

Note that in the usage of equation (168]) there is another limitation of close encounter, which 
has already been addressed in Section 13.31 



5. Discussion: Dust Gap Opening 

5.1. Order-of-Magnitude Estimate of Dust Gap Opening Criterion 

Paardekooper and Mellema (2004, 2006) showed that dust gap may be opened up even 
if the planet mass is too small to open up a gap in a gas disk using numerical calculation. 
We reconsider this result from the analytical point of view in this section. 

Since there is no long-term radial motion of fluid elements, it is essential in dust gap 
formation that gas and dust moves in a different velocity. Terms that are important in 
equation (168]) are therefore those proportional to z/f2 p /(z/ 2 + Q^) for small particles. 

We first consider the disk without global pressure gradient. In this case, the motion 
of small dust particles are mainly dominated by the third term of equation (168]) when they 
are as close to the planet as b <^ 2H and by the last term when b is larger. The third term 
is the attraction of particles by the planet's gravity and they migrate towards the planet. 
The last term on the other hand originates from the axisymmetric flow structure around the 
planet. This term causes dust particles migrate away from the planet. In either case, dust 
particles around the planet with b ~ H forms a gap around the planet. Taking b ~ H, very 
rough order-of-magnitude estimate (neglecting numerical coefficients of order unity) gives 
the migration timescale of dust particles as 

where Tk denotes the Kepler timescale and M c is the mass of the central star. This is 
the timescale of dust particles changing its orbit by the order of scale height and gives the 
timescale of the formation of dust gap with width of the order of the scale height. 

If there is global pressure gradient, the first term of equation (168]) may dominate over 
the terms considered above. Since global pressure gradient causes systematic inward or 
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outward motion of particles, dust gap around the planet is not expected to form if this term 
dominates. The timescale of orbital change of dust particles by the order of scale height due 
to the first term of equation (|68|) is 

If the magnitude of rj is of order (H/r p ) 2 ~ 10~ 3 , which is a standard value, dust gap may 
not appear around the planet since r v < r gap . However if rj is as small as 10~ 5 , dust gap 
formation may be possible. The value of rj depends on disk model and can take values of 
wide range. From theoretical point of view, small value of r\ is preferable in order for meter- 
size particles to survive in a Minimum Mass Solar Nebula. One mechanism of decreasing 77 
around ice line is suggested by Kretke and Lin (2007). 

Comparing the first and the third term of equation (|68|) , we obtain the criterion of dust 
gap formation. Again, for b ~ H, and neglecting numerical coefficients of order unity, we 
have, for the gap formation criterion, 

1 * 1. (90) 
M p r p 

We note that this criterion does not depend on drag coefficient of the particles. However, 
the implicit assumption is that gravitational scattering, which is the fourth term of equation 
( |68l) . is ineffective, which is proportianl to Vt^/{v 2 + Q 2 ). Therefore, this condition applies 
to particles with v ^ fl p . Taking rj ~ (H/r p ) 2 and H/r p ~ 0.05, this gives M p /M c ^ 10~ 4 , 
which explains the condition obtained by Paardekooper and Mellema (2004, 2006). 

Timescale constraint of dust gap formation may be obtained by comparing the above ob- 
tained timescale with the timescale of planet migration or gas dispersal timescale. Timescale 
of the type I planet migration of the radial distance comparable to the scale height is given 
by Tanaka et al. (2002) as 

~ a »w^ U • (91) 

Therefore, comparing this with r gap , we have 



r P Erg v 2 + Ql 
H M c vVt p 



& I- (92) 



Dust gap formation is efficient at late stage of planet formation when gas starts to dissipate 
and type I migration timescale becomes long. If H/r p = 0.05 and Hr 2 /M c = 10~ 3 , this 
criterion gives v/VL p <^ 50. If there is no type I migration, on the other hand, this constraint 
is relaxed and dust gap formation constraint is given by the timescale of disk gas dispersal, 
r gas . If r gas is of the order of 10 6 years, gap of particles with v /VL p <; 10 4 may be opened up. 
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5.2. Model of Radial Distribution of Dust Particles 

In order to investigate the qualitative behavior of dust particles around the planet in long 
time, we model the motion of particles in the disk by one-dimensional advection equation 

^%^ + Wb [Vbm{t > b)] = ° (93) 

where Vb{b) is the radial velocity of the dust particles whose semi- major axis difference 
with the planet is b and N(t, b) is the appropriately normalized number of the particles 
with semi-major axis difference b at time t. We use equation (I68I) for and solve (1931) for 
10 6 yr for particles with various drag coefficients. We have used second-order scheme with 
monotonicity condition (van Leer 1977). We fix H and c for simplicity in our computational 
domain and a planet with GM p /Hc 2 = 10 _1 is fixed at the origin. Planet mass corresponds 
to 2M e for H = 0.05AU and c = 10 5 cm/s. The computational domain is r# < b < 10H and 
dust particles are homogeneously distributed initially: Nq = N(t — 0, x) — 1 in this region. 
We have assumed that there is no dust particle out of this domain. For T that appears in 
equation fl68l) . we have assumed 

T = 27r/\Q K (r p + b)-Q K (r p )\ (94) 

where Qk is Keplerian angular velocity and r p is the location of the planet. 

Figures [T7] shows the snapshots of the distribution of dust particles of various drag 
coefficients at tQ p = 5 x 10 3 , tQ p = 10 5 , and tQ p = 10 6 for 77 = 0. The gap of particles with 
v ~ lOfip has been opened up after several 10 5 Kepler time. The gap of particles v <^ Q p may 
be the artefact of boundary condition, since the radial velocity of particles at the nearest 
edge of the planet is positive but we have assumed there is no dust particles outside the 
calculation domain. However, there is no boundary effect for small particles with large drag 
coefficient, since radial velocity of dust particle is negative at the inner edge and positive 
at the outer edge. We note that since it takes about 10 5 years for dust gap formation for 
the planet with 2M e , this effect is small if the gas pressure gradient is significantly large, 
i] > 0. However, for larger mass planets, planet perturbation on the dust motion is large 
and dust gap formation is possible even in the presence of pressure gradient, as seen in the 
order-of-magnitude estimate provided in the previous section. 

Finally, we briefly discuss the observational implication of dust gap in the disk with 
zero pressure gradient. The width of the dust gap that forms around a low mass planet is 
several times the scale height. Therefore, it might be challenging to detect a low mass planet 
embedded in a gas disk at 1AU. However, if there is a planet at several tens of AU, the 
width of the gap can be of the order of several AU, which can be resolved in near future. 
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The gap formation timescale in this case is of order 10 5-6 years, which might be comparable 
to disk life time. Therefore, it may be possible to find a low mass planet by observing a dust 
gap in a gaseous nebula, even if the planet itself is too small to be observable. 

6. Summary and Future Prospects 

In this paper, we have investigated the motion of particles embedded in a gas disk in the 
presence of a small mass planet. This is the full analytic calculation of the particle motion 
including the gas disk and the planet considering the non-axisymmetric pattern of the gas 
around the planet produced by the gravitational interaction between the planet and the gas. 

Our main result is equation (168|) which describes the change of the particle semi-major 
axis, which has been in circular motion initially, after one distant encounter with the planet. 
In addition to well-known particle migration towards the central star due to the pressure 
gradient (the first term) and the gravitational scattering by the planet (fourth term), we have 
derived (1) the effect of steady mass accretion (the second term), (2) attraction of particles 
towards the planet due to the planet's gravity (the third term), and (3) contribution of the 
structure of the gas disk produced by the planet's gravity (the fifth term). We find that 
only the axisymmetric structure of the gas contributes to the semi-major axis evolution of 
particles and the spiral structure is ineffective in the absence of vorticity source. 

All the terms considered, in the absence of the global pressure gradient and steady 
accretion flow towards the central star, we find that (see Figure [3]) (1) particles with v ^ Q p 
will migrate towards the planet if they are close to the planet because of the gravitational 
attraction by the planet, but they migrate away from the planet if they are far away, primarily 
due to the effect of the gas structure, (2) large particles with v <C Q p will be scattered away 
from the planet because of the gravitational scattering, and (3) there is a parameter range 
where particles are scattered away from the planet when they are close to the planet while 
attracted towards the planet if they are located at far away, thereby particles are accumulated 
at an equilibrium distance from the planet (i.e., shepherding by gas drag and gravitational 
scattering). This shepherding mechanism acts for particles in both interior and exterior 
orbits as long as the radial pressure gradient of the gas is negligible. 

We have checked the validity of our formula by solving three-body problem with and 
without Hill's approximations. Our main assumptions in deriving equation (1^5]) are (1) 
initially circular orbit, (2) distant encounter, and (3) local approximation. Validity of the 
first assumption is justified when the eccentricity of the particle is damped within the time 
taken for each synodic encounter that is given by equation fl51|) . The valid region for the 
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second assumption is shown by the contour of small deflection in Figure [TTJ The condition 
for the third assumption is simply given by equation flH7|) . Hill's approximation is good as 
long as we consider the orbits in the vicinity of the planet. In short, our formula is especially 
useful in predicting the motions of dust particles with stopping time t e Q p <; 10 2 . For particles 
with smaller drag, our treatment predicts only quantitative behavior of non-resonant orbital 
evolution of particles whose semi-major axis is close to that of the planet. We also note that 
since we have made use of linear calculation for the flow, our treatment is restricted to low 
mass planets with the mass up to the orders of several tens of Earth mass. 

Using equation fl68|) . we have discussed the criterion for dust gap formation. The con- 
dition depends on the value of global pressure gradient, mass of the planet, and disk scale 
height. It is possible to derive the condition suggested by Paardekooper and Mellema (2004, 
2006), and we have derived more general criteria. We have calculated the qualitative long- 
term behavior of dust particles around a planet with 2M e . It is indicated that the gap 
of small dust particles with width of the order of several scale height may be opened up 
after several 10 5 Kepler time. This long timescale has not been reached by full numerical 
simulations yet. 

We have simplified the description of the drag force by assuming the force proportional 
to relative velocity with the gas. For large particles, however, the drag force is proportional to 
the square of the relative velocity, which is not treated in our calculation. We have also made 
a simplification by considering two-dimensional problem. In considering the vertical motion 
of particles, it is necessary to do three-dimensional calculations. The effects of magnetic field, 
self gravity, and turbulence of the disk will be the subjects in this line of analytical study. 
It is also of interest to extend our formulation to include the finite box size and investigate 
how the resonance effect is described in Hill's system. 
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A. Derivation of Equation ( 120 



In this section, we briefly show the derivation of (120]) . The method is the same as that of 
Goldreich and Tremaine (1980), Henon and Petit (1986), or Hasegawa & Nakazawa (1990). 

We consider the restricted three-body problem. In other words, we consider the problem 
with v = 0. In this case, Jacobi energy Ej conserves throughout the particle orbit. Jacobi 
energy is given by 



Ej = -r 2 p nl(h 2 + k 2 ) - + (Al) 



1 

2 

where ip p is the gravitational potential of the planet. Therefore, orbital elements before and 
after the encounter with the planet are related by 

- l^b 2 = ^(Moo) 2 + Moo) 2 ) - 3 -n 2 p (b + Ab) 2 . (A2) 



Assuming Ab -C b, we obtain 



2r 2 

Ab ~ (h(oo) 2 + k(oo) 2 ) . (A3) 



On the other hand, from equation (JUJ), A6 is given by 

c\ POO 

i4 P J -oo 



Therefore, we obtain 

r 2 fl 



/oo r ^() 
F y {t)dt = ^(h(oo) 2 + k(oo) 2 ). (A5) 
36 

From equation ( fT2l) and ( [TBI , the elements h and k after the encounter are given by 

1 r°° 

h(oo) = — / du {F x (u) sm[Q p u] + 2E y (u) cos[Q p u}} (A6) 

r p"p J -oo 

1 f°° 

k(oo) = — — / du{F x (u) cos[Q p u] — 2F y (u) sm[Q p u]} . (A7) 

r p"p J-OO 

Approximating the trajectory of the particle by circular orbit, we have 

FM) ~ -sgn(6)^^ — (A8) 

& 2 (l + (9/4)(fi p t)2) 3 / 2 

/n GM p (3/2)0,* , 4 „, 

F y (t) ~ sgn (&) — * v 1 1 p rjr. A9 

& 2 1+ 9/4 O p i 2 3/2 
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Integrating (IA6I) and (IA7I) . we have 



fc(oo) = — sgn(fe) 



/i(oo) = 
8 1 GM, 



9r p Q2 &2 



#1 



+ 2K 



Substituting these results into equation ( 1A5I) . we finally obtain 



2/\/f2 r 



F y (t)dt 



64 g 2 M^ 
243 6 5 fi3 



2^n 



(A10) 
(All) 

(A12) 



We note that direct substitution of equation (1A9|) into equation (IA4I) results in zero, indi- 
cating that higher order of the expansion is essential. 



B. Proof of Equation ( 164 



In this section, we show that the right hand side of equation fl63l) vanishes up to second 
order perturbation and prove equation fl64|) . The equation we prove is 



d r T dSv (1) r h 

n dt5v ( i ) —5v^ ) + J dh-^-ixcycih^j dt 2 6vi 1 \x e ,y e (t 2 )) 



-I dh 
/o oy 



x C} y c {ti)) \ dt 2 5v v 1) (x c ,y c (t 2 )) = 



We Fourier transform the perturbation in the ^-direction, 

5f(x,y)=Jf+J2 5 ~fk y (x)j k 

ky^Q 



(Bl) 



;b2) 



where 5f denotes any perturbation quantity, 5f is the average over y, and f ky is the Fourier 
component of non-axisymmetric modes. Since physical quantities are real, 5f ky = of_ k . In 
this section, we drop superscript (1) since all the perturbed values refer to first order values. 

Since we integrate over the circular orbit, integration with respect to t is converted to 
that with respect to y by 



y 



3 

2 p 



(B3) 



It is possible to show that the right hand side of equation ( IBlj) is equal to 



8L, 



.Re 



$ v yk 

dSv 



dx 



8L y y ii m 

9^2 l^k y >o ky ™ 



dSv xk 



dx 



y -5v 



8L, 



ky V ky 



+ 



4L„ ^_ 



5v y 5v x (x c ,L y /2). 



(B4) 
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The last term is zero as long as the box size is large enough. Using equation of continuity 

3 d 
- -ikyxQp^-^ + -j^.S v xk y + iky8v y = (B5) 

and conservation of vorticity 

' - ~ 1 5S fc 

ik y Sv xky - -fip^^ = 0, (B6) 



—5 v. 



dx y ^ 2 " >;„ 

we can show that the first three terms cancel. Hence, equation is proved. 



C. Mass Flux of Sound Wave Propagating in a Homogeneous, Static Medium 



In this section, we demonstrate that the sound wave that propagates in a homogeneous, 
non-rotating medium carries mass flux, and therefore, show that the vanishing mass flux 
given by equation (|65|) is not a general conclusion but specific to the problem in consideration. 
We consider one-dimensional system which extends in x > 0, and at x — 0, there is a forcing 
that creates a sound wave, which is turned on at t — 0. The full system of equations are the 
equation of continuity and Euler equation in one-dimensional system 



i + a~J px) = 

dv dv (? do 

— + v— = — jf- + Scos{nt)6 D {x)e{t), 

at ox p ox 



(CI) 
(C2) 



where p is density, v is velocity, c is sound speed, S is the amplitude of the source, Q is the 
frequency of the forcing, 8d(x) is the Dirac's delta function, and 6{t) is step function. We 
consider background state with constant density p and vanishing velocity. 



The first order fluctuation caused by the forcing is given by 

^7 + 7T 6v = 
at po ox 

+ c 2_^_^ = S cos(ttt)5(x)6(t) 
at ox po 

and the second order perturbation is caused by the first order perturbation as follows 

V 1 ) 



at po ox ox |_ Po 



Ufa® + = — 

dt dx po 2 dx 



8pV> 
Po 



8vW 

2 



(C3) 
(C4) 

(C5) 
(C6) 
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Mass flux Fm = pv is given by, up to the second order, 

F M = p Q 5v {l) + Po 5v^ + 5p {1) 5v {1) 



(C7) 



Since the sound wave propagates to x — > oo and there is no reflection, the appropriate 
solution for the first order is given by 



6pW 
Po 

6vW 



S 



s 



cos(flt — kx), 
cos(f2t — kx) 



where k = Q/c and the boundary condition at x = is considered 
Equation for second order velocity fluctuation is given by 

d 2 



c 2 dt 2 dx 2 



Sv (2) = J&P cos [2(nt-kx)}. 



Since all the perturbation must vanish at t — 0, we obtain 



5v {2) = -—fit sin \2(nt - kx)} . 



Therefore, mass flux is given by 
Fm 



(C8) 
(C9) 



(CIO) 



(Cll) 



p s 2 + poS_ cog ^ _ ^ + sin j 2(fit _ kx ^ + cos p( nt _ kx y _ (C12) 



2c 3 



2c 3 



2c 3 



Therefore, there is a positive and finite mass flux p S 2 /2c 3 remains even after spatial and 
temporal averages are taken. It is also possible to show that a fluid element move to x > 
on average by calculating path line. 



D. Interpolation Methods 

In the numerical calculation of large drag coefficient, stopping time of the particle u^ 1 
is very small that the particle remains in the same grid cell of the hydrodynamic calculation 
for several time steps. In this situation, it is necessary to have a smooth data of gas velocity 
even in sub-grid scale in order to obtain a smooth results. In this section, we describe how 
this can be realized. 

Our method of interpolation is to consider the value of one grid cell as a representative 
value of physical data around the grid cell with a certain width, and add the contribution 
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from nearby grid cells to find the value of physical quantities at a location of interest. Let 
f(x) be the value of a physical quantity of interest and the data of f(x) is given at discrete 
set of grid cells at x — Xj. We denote the data of f(x) at x = xi by = /(x,) and A by 
grid size A = x i+ i — x,. Let the grid cell which is closest to the location of interest x be Xio. 
Our method of interpolation approximate the value of / at x by 



/(*) ~ "T^oTT — ( D1 ) 



where J is an integer and 77 is a numerical factor. If rj is large, we have a more smoothed but 
more damped value. We find I = 20 and 77 = 4 gives a reasonably smooth results. 
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dust size (cm) 

Fig. 1. — An example of the relationship between the dust size and the drag coefficient v used 
in our numerical calculations, given by equation ([6]). Reciprocal of stopping time, v = t~ x 
is shown as a function of dust size in centimeters for gas densities with p = 10~ 8 g/cm 3 , 
p = 10~ 10 g/cm 3 , and p = 10" 12 g/cm 3 . 
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Fig. 2. — Comparison between the azimuthally averaged profile of L y 5v y /c obtained by the 
numerical calculation (solid line) and analytic expression floTl) (dotted line). The results of 
the numerical calculation is normalized by GM P /Hc 2 . 
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Fig. 3. — Contour plot of Ab/T given by equation ( 1681) normalized by iff2 p for a planet 
with GM p /Hc 2 = 10 -2 . The horizontal axis shows the distance from the planet and the 
vertical axis shows the drag coefficient of the particle. Global pressure gradient and steady 
mass accretion is neglected (77 = ( = 0). We note that our analytic approach is limited 
by assumptions of local approximation, initially circular orbit of the particle, and distant 
encounter. This figure has only qualitative meanings for particles with u/Q p <^ 10~ 2 , with 
initial semi-major axis difference b <^ 3rn for u/Q p ^ 1 — 10, or with b <^ th for u/Q p ^ 10. 
Detailed discussions in the limitation of analytic calculations are given in Sections 13.31 and 
IU and are summarized in Section [6j 
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Fig. 4. — Evolution of b for a particle with zero drag coefficient initially located at b = H. 
Symbols show the values of b obtained at the box boundary. The value of b during the 
encounter is also indicated by dashed line upto t = 250f2~ 1 . 
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Fig. 5. — Evolution of b for a particle with zero drag coefficient initially located at b- m [ = 3 AH 
(plus) that corresponds to I = 5 case in equation (ITT)) and &i n i = 3.8H (cross) that does not 
satisfy equation (jTTj) . Vertical axis indicates the value of (b(t) — b ini )/H. The value of b(t) is 
obtained when the particle reaches the box boundary. Note that vertical axis of this figure is 
in logarithmic scale and the values smaller than 10 -9 are not plotted. Horizontal axis shows 
time. 
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Fig. 6. — Evolution of b for particles initially located at x = H with different drag coefficient. 
Planet mass corresponds to GM P /Hc 2 = 10~ 2 . Particles with small drag coefficients shows 
oscillation, while particles with large drag coefficient shows a systematic decrease or increase 
of semi-major axis. 
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Fig. 7. — Comparison between the local numerical calculation (solid line) and analytical 
formula (dashed line) of Ab/T for particles with u/Q p = 10~ 2 . The value of Ab/T is the 
average of a number of encounters and calculated using equation (1791) . The gravitational 
force by the planet and the spiral pattern of the gas around the planet are both considered. 




Fig. 8. — Same as Figure [7] but for particles with u/Q p = 1. 
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Fig. 9. — Same as Figure [7] but for particles with u/Q p = 10 2 . 
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Fig. 10. — Long-term evolution of the semi-major axis of the particles with v = Q p . The 
mass of the planet corresponds to GM p /Hc 2 = 10~ 2 . The calculations are done for the disk 
model with zero pressure gradient and zero accretion flow, 77 = ( = 0. 
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Fig. 11. — Deflection of the particle's semi-major axis after the first encounter with the 
planet. The amount of deflection is quantified by equation (1551) . Horizontal axis denotes the 
particle's initial semi-major axis and vertical axis denotes the particle's drag, or equivalently, 
size. Initial semi-major axes are shown in terms of scale height (bottom axis) and Hill's radius 
(top axis). Gray scale shows the amount of deflection, and the contours of 10% (solid line) 
and 1% (dashed line) deflection are shown. Particles in black regions are trapped by the 
planet (that is, mutual distance between the particle and the planet becomes smaller than 
half of Hill's radius) or in horseshoe orbit. 
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Fig. 12. — Examples of the orbit of the particle when they encounter with the planet. 
Particles with is/Q p = 10~ 4 (solid line) and v/Vt p = 10 3 (dashed line) are shown. Both 
particles are initially located at b = 2rn- Only the region (lrn < x < 2rn, — 50rn < 
y < 50r H ) is shown. The particle with small drag is strongly perturbed by the planet and 
eventually captured within the Hill's radius, while that with large drag can escape owing to 
the drag by the background gas. 
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Fig. 13. — Comparison between the global numerical calculation (solid line) and local calcu- 
lation (dashed line) of Ab/T for particles with v/£l p = 1CT 2 . The average of the change of the 
semi-major axis of the particle after 100 encounters with the planet is shown. Horizontal axis 
shows the initial location of the particle normalized by disk scale height, b/H = (r — r p )/H. 
We use disk aspect ratio H/r p = 0.04 and the mass ratio between the planet and the central 
star is M p /M* = 6 x 10~ 7 . The disk aspect ratio is chosen in such a way that the value of 
L y we have adopted in the integration of Hill's equation corresponds to 2nr p , and the planet 
mass is chosen in such a way that GM P /Hc 2 = 10~ 2 , which is also the value we have used 
in local calculations. Keplerian rotation of the gas is assumed so i] = ( = and there is 
no modification of the gas structure by the planet's gravity. Note that global calculations 
and local calculations are in good agreement except for the location of the resonances. The 
position of 4 : 5, 5 : 6, and 6 : 7 resonances are indicated by arrows. 
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Fig. 14. — Same as Figure fT3l but for particles with v/Q p = 1. Particles initially located at 
b/H < 1.4 has been trapped in the vicinity of the resonance so the average rate of orbital 
change over 100 encounters is small. See also Figures [10] and [16] for the time evolution of 
semi-major axis difference b for particles with b/H <; 1.4. 




Fig. 15. — Same as Figure fT3l but for particles with u/flp = 10 2 . 
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Fig. 16. — Evolution of semi-major axis of particles with u/Q p = 1 in global run with 
7] = ( = and the gas is assumed to rotate at Kepler velocity. Osculating elements of 
the particle at the opposition point is plotted. Particles in both interior and exterior orbits 
migrate towards the planet and trapped in the orbit of ~ 0.3H. This figure corresponds to 
Figure [10] that shows the results of local calculations. 
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Fig. 17. — Time evolution of the distribution of dust particles around 2M® embedded in 
a solar nebula at tfl p = 5 x 10 3 (top left), tQ p = 10 5 (top right), and tQ p = 10 6 (bottom 
left) with various v. The horizontal axis denotes the semi-major axis difference between the 
planet and the particle. Gray scale indicates the fraction of the remaining particles. The 
pressure gradient factor r\ is set to be zero. 



